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Abstract 

The surface or contact potential at the water Hquid-vapor interface is discussed in relation to 
determinations of absolute ion hydration free energies and distributions of ions near the interface. 
It is shown that, rather than the surface potential itself, the net electrostatic potential at the center 
of an uncharged solute can aid both in relating differences between tabulations of hydration free 
^0 ■ energies and in explaining differing classical and quantum surface potential estimates. Quantum 

mechanical results are consistent with conclusions from classical simulations that there is a net 
driving force that enhances anion density at the surface relative to cations. 
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I. INTRODUCTION 



Both the sign and the magnitude of the surface potential (psp across the water hquid/vapor 
interface have eluded experimental and theoretical consensus, although recent progress has 
been made. Electron holography experiments suggest a surface potential of +3.5 V for 
vitrified ice Two density functional theory (DFT) studies have obtained values of +3.1 
V ^ and +3.6 V ^ for the water liquid-vapor interface, consistent with the experimental 
result. A related DFT study revealed a shift in ion HOMO states in clusters compared with 
bulk periodic boundary calculations, suggesting a surface potential value of +4 V [3]. Such 
a large potential change across a relatively narrow interface would naively suggest enormous 
electric fields that would have substantial chemical effects. Surface potentials computed with 
classical point-charge or polarizable models, on the other hand, are typically in the range 
-0.6 to -0.5 V liij. (The quoted values are for the SPC/E, TIP3P, TIP4P, Dang-Chang 
(polarizable) and TIP4P-FQ (polarizable) water models; recent work has shown that the 
surface potential of the TIP5P water model is closer to —0.1 V [l^.) 

It is generally accepted that the surface potential is not directly accessible to thermo- 
dynamic measurement [ll|. Nonetheless, electrochemical measurements have been used to 
indirectly infer surface potential values with magnitudes much smaller than 3.5 V, 



ranging from —0.4 V to +0.3 V, as summarized by Randies [18|. The experiments have 
employed a variety of approaches, including models that involve Born theory to estimate 
bulk hydration contributions for large ions [l^. Those estimated bulk contributions are 
subtracted from the so-called real hydration free energies (the free energy terminology is dis- 



cussed below), allowing estimation of a surface potential contribution (—0.3 V in Ref. [19 
Alternatively, Randies discusses studies of the temperature dependence of cell emfs [jj 
early experiments suggested a small positive surface potential value (0.1 to 0.2 V) j20|, but 
the final result quoted in Ref. [l8| is 0.08 ± 0.06 V, or effectively zero. Subtleties related to 
the effects of interfacial potentials on single-ion thermodynamic properties are discussed in 
the recent books of Fawcett j2l| and Hiinenberger and Reif [2^ . 

How can the widely differing electron holography and electrochemical results be related? 
Kathmann et al. jlj have taken an important step toward answering this question by showing 
that, with omission of certain regions of space (namely the space inside the water molecules) 
in computing average potentials, the computed and electrochemically estimated surface po- 
tentials come into closer agreement (as do the classical and quantum mechanically calculated 
results). In related work involving classical point-charge and polarizable simulations. Harder 



al p 

and Roux [8| and Vorobyov and Allen 2^ have shown that it is the real hydration free 
energy, as opposed to the intrinsic free energy (also defined below), that provides relatively 
consistent results due to a cancellation of model-dependent molecular quadrupole contri- 
butions between the solute-solvent and water liquid-vapor boundaries. Two goals of the 
present paper are to extend the analysis in Refs. [8| and |23| to quantum calculations and to 
provide a firmer statistical mechanical resolution of the above question. In the process, the 



surface potential as defined by Landau, Lifshitz, and Pitaevskii [2J], and that inferred from 
electrochemical experiments [l8|, are shown to be different quantities. 

The discussion presented here impacts on both measurements of single-ion absolute hy- 
dration free energies and the distributions of ions near the water surface. It will be shown 
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that the proton (real) hydration free energy estimate obtained by combining bulk and cluster 
data, and using the cluster-pair approximation 25|] , does include a full contribution from the 
bulk water surface potential under the assumption that the common intersection behavior 
observed in the data analysis continues to hold for clusters larger than n = 6. Differences be- 
tween the proton hydration free energy reported in Ref. 25|] and other tabulations involving 
extra-thermodynamic assumptions 26|, |27(] are largely due to omission of the electrochemical 
surface potential, which is shown below to be closely related to the net potential at the 
center of a neutral, ion-sized particle. For typical monatomic ions, that net potential shows 



limited size dependence |28H31 



and persists up to relatively large solute sizes |29| . 
In previous classical studies, it has been shown that the net potential influences the 
distributions of ions near the water surface [ol, [32| . But the question remains whether this 
influence persists when water is modeled at the quantum mechanical level [s^. The statistical 
mechanical basis for addressing this question is first developed, and results are then presented 
for an ion in a water droplet (of sizes n = 105 and n = 242). Again, the net potential at 
the center of a neutral particle is shown to play a prominent role. The quantum calculations 
are performed on clusters at the all-electron DFT and Hartree-Fock levels, thus eliminating 
potential complications arising from pseudo-potentials and intricacies of the Ewald potential 
in periodic boundaries. The classical droplet results (for n = 105 and n = 511 clusters) are 
similar to the previous classical simulations in the slab geometry [9|. The quantum results 
for the net potential are qualitatively similar to the classical point-charge (SPC/E water 
model) results. As in the classical studies discussed above, the reason for the similarity of 
the net potentials (in spite of the enormous difference between the computed water surface 
potentials) is due to the cancellation of a quadrupole contribution that differs widely (in sign 
and magnitude) between point-charge and quantum water models. In the quantum case the 



canceling terms are significantly larger. 



II. OUTLINE 



The main focus in this paper will be on the simple situation of a single ion inside a water 
droplet. The bulk (infinite dilution) limit is attained as n — )■ oo, where n is the number of 
water molecules. To make contact with previous studies and recent experiments js^ (on 
a related large, chaotropic ion, SCN~), the I~ ion is primarily examined in the calculations. 
The hydration free energy of the Na"*" ion is also computed. It will be assumed that the ions 
are monatomic and monovalent to ease the notation. 

The discussion begins with a review of the electrostatics of the surface potential and its 
decomposition into molecular dipole and quadrupole contributions. Then two forms for the 
ion hydration free energy are discussed that will be useful in analyzing surface potential ef- 
fects. Next, statistical mechanical expressions are derived for the hydration free energy, and 
connection is made to a simple and robust approach for extracting bulk single-ion quantities 
from ion hydration data. Two exact approaches for partitioning the free energy are summa- 
rized; the partitioning aids both in relating to experimental measurements and in isolating 
the role of the surface potential in the free energy. 

The discussion then turns to recent experimental determinations of the absolute hydration 
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free energy of the proton [25|] , and adsorption enthalpy and entropy changes for a large anion 
approaching the water liquid- vapor surface j33|. At first glance, these experiments would 
appear to be unrelated. It is shown, however, that the same free energy shift included in the 
proton hydration free energy measurement leads to a net driving force that enhances anion 
density relative to cations at the water surface. Results of classical and quantum calculations 
on ions in water droplets are presented that provide preliminary support to this view. 



III. THE SURFACE POTENTIAL 

The potential through a liquid-vapor interface (starting in the vapor, and passing to 



the liquid, I) can be obtained by integrating the one-dimensional Poisson equation (34 



(j){z) - <j){v) = 47r / {z' - z)p,{z')dz' (1) 

J V 

where z is directed along the surface normal, and Pq{z) is the total average charge density 
at a particular z location. In the macroscopic limit, integrating deep into the liquid, the 



resulting surface potential is [24 



■'sp 



4:71 / zpq{z)dz (2) 



Thus, the surface potential jump is the result of a dipole layer. That surface dipole in turn 
results from the average distribution of all the charges, however, and does not simply reflect 
molecular dipole orientations at the interface. In simulations, it has been shown that the 
potential shift across the water surface occurs on a length scale of less than 10 A jl, 0]- 
In pioneering work, Wilson, Pohorille, and Pratt 35, 3^ showed that the surface potential 



of water computed from the molecular dipole contributions alone yields incorrect results. 



Building on the development of Jackson [37[ for the average macroscopic charge distribution 
in terms of the molecular multipole contributions, they showed that there is also a substantial 
contribution to the surface potential from the molecular quadrupoles: 



(t)sp = 4tx I P,{z')dz' - 47r [Q,,{zi) - Qzz{zv)] (3) 

where Pz is the z component of the molecular dipole density, and Qzz is the zz component 
of the density of molecular quadrupoles. The separate dipole and quadrupole contributions 
depend on the choice of molecular center, but their sum does not. The multipole expansion 
truncates sharply at the quadrupole term since it is a gradient expansion, and the two 
integrations leave spatial derivatives for the higher-order terms that are zero in the bulk 
phases. 

The density of quadrupoles in the vapor is negligible, and the quadrupole density in the 
liquid is given by the trace of the quadrupole tensor: 

Qzzizi) = ^Tr [Q^,{zi)] = (^Y. ^i^i - ^m) E ) (4) 
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for a molecular point-charge model [6|. The m label is for the molecules, while i labels the 
charges within a molecule. For the quantum mechanical case, the charges are replaced by 
the continuous charge density, and the sum over i is replaced by an integral; the general 
formula for the quadrupole matrix elements is then [37| : 

Qa^ = sJ Xy^p{^')£x' (5) 



where Xa are the Cartesian x, y, z components and p(x) is the charge density. Jackson [3j 
provides a careful discussion of why the nonzero-trace form of the quadrupole is required 
for the expansion of the charge density, as opposed to the traceless form (the trace adds a 
required / = term to the multipole expansion). 

An interesting feature of these results is that the quadrupole contribution is a property 
of the bulk liquid and has nothing to do with molecular orientations at an interface. We will 
see below that the classical point-charge quadrupole trace for water is positive, leading to a 
negative contribution to the surface potential, while the quantum result is large and negative, 
producing a large positive contribution. Another point is that, if we consider the net potential 
inside an ion-sized cavity embedded deep in a water droplet, the quadrupole contribution 
cancels since the quadrupole densities both in the vapor and inside the cavity are zero. The 
process of moving a non-interacting test charge from the vapor, into the droplet, and then 
into the cavity results in the crossing of two 'interfaces', the water liquid- vapor interface, and 
the water-cavity boundary, leading to the cancellation. Then the net potential at a cavity 
center reflects a dipole contribution. But that dipole contribution arises from differences in 
water orientations between the solute-water boundary and the water liquid-vapor interface. 



IV. ION HYDRATION FREE ENERGY DEFINITIONS 



The above discussion views the surface potential as an electrostatic effect involving the 
average charge distribution through the interface. The equilibrium electrochemical potential 
fi for a given ion is typically written as [ll|, [sF 
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■'sp 



(6) 



where k is Boltzmann's constant, T is the temperature, p is the ion number density, A is 
the thermal de Broglie wavelength, /i*^^ is the real hydration free energy, pf^^ is the intrinsic 
hydration free energy (in the absence of an interface) j4o|, and q is the ion charge. The 
real hydration free energy p'^^ is the free energy change for moving an ion from a vacuum, 
across the water surface, and into the bulk. The surface potential is the difference of the 
inner (Galvani) and outer (Volta) potentials j38|. In what follows we will assume that the 
potential is zeroed in one of the phases, here taken as the vapor phase. Both /i*^^ and (psp can 



■'sp 



-80.7 kcal/mol 



be obtained experimentally (/i*^^ = —59.3 kcal/mol for I~ 2^ while 

estimated from electron holography [l|), resulting in p^^^ = -1-21.4 kcal/mol. These results 
show that pf^j differs from any electrochemical estimates of bulk hydration free energies 
[sl. Ell. (The DFT calculations of ion hydration free energies in Refs. 113 included a 
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procedure to subtract off the large quadrupole contribution arising in periodic boundaries; 
that procedure leads to free energies more in line with the bulk values discussed below.) 

To set the stage for the discussion below, and following Harder and Roux [1], we choose 
to re-define the real hydration free energy as 

Ai'" = A^r + (l(t>np (7) 

where /i^^ is a bulk hydration free energy, while 0„p is the net potential at the center of an 
uncharged solute. The rearrangement amounts to adding and subtracting a term which is 
the charge times the local potential due to nearby waters sampled with an uncharged solute. 
As discussed above, this leads to a cancellation of quadrupole terms. With the results 
obtained below, the resulting I~ value is —70.9 kcal/mol (using the value (pnp = —11.6 
kcal/mol-e or —0.50 V). With this re-definition, is much closer to estimates based on the 



parameterization of Latimer, Pitzer, and Slansky [43| (/x^^ = —71.4 kcal/mol in Ref. [44 
see below). We note here that all hydration free energies listed in the present paper refer to 
a transition from a 1 M vapor to a 1 M solution phase concentration. Thus, for example, all 
free energies taken from Refs. j25| and [3] are uniformly shifted downward by 1.9 kcal/mol. 

It is clear that surface potential estimates indirectly obtained from electrochemical mea- 
surements (and perhaps further assumptions) are close in magnitude to the (pnp value listed 
above (—0.50 V). The cost of the rearrangement, however, is a slight solute-size specificity to 
(/)„p (of ma gnit ude 1 kcal/mol-e variation for neutral particles sampled with Lennard- Jones 



potentials) [28|, l30|, |31 



To summarize, six quantities have been introduced: the electrochemical potential /i, the 
real hydration free energy /x^^, the intrinsic hydration free energy /i^^^, the bulk hydration 
free energy /i^^, the surface potential (pgpi and the net potential at the center of an uncharged 
solute (j)np- The terminology used here was chosen to follow the existing electrochemical and 



simulation literature |38|, |40| . All of the excess quantities and the potentials can be estimated 



from experiment, but the surface potential cannot be measured thermodynamically. 



V. STATISTICAL MECHANICS OF ION HYDRATION 



We take the Potential Distribution Theorem (PDT) as the fundamental expression for 



the chemical potential [45|, |46 



/i = fcTln pA 



where 



^ex ^ (exp(-e/A;T))o 



(9) 



and e is the interaction energy of the ion with the surrounding medium. The '0' subscript 
on the thermal average implies the sampling is conducted with no ion present. Eq. M is 
equivalent to the inverse form 



/x"^ = kT\n {exp{e/kT)) 



(10) 
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where now the samphng includes the ion. The excess chemical potential in Eqs. M and [TU] 
includes any possible contribution from the surface potential since that contribution arises 
from the interaction energy of the ion with a distant dipole layer. 

Above, the uncoupled sampling in Eq. [Hlhas no ion present; this is a correct formula for 
the excess chemical potential, but can lead to computational difficulties for larger particles 
due to the hard ion core. Thus in what follows we will assume that the free energy has 
already been estimated for the uncharged ion, and the '0' subscript will imply sampling with 
the uncharged ion present but no electrostatic coupling. The 'es' subscript indicates the 
electrostatic contribution to the energy or free energy. Then the total hydration free energy 
is /i^^ = fi^dw + f^esy the sum of the non-polar and electrostatic contributions. 

The resulting electrostatic contributions to Eqs. M and [TU] can be rewritten exactly as 

^iTs = {ees)o - kT\n (exp [-(^e. - (ee.)o)] /kT))^ (11) 

and 

fit: = i^es) + kT\n (exp [(e^s - (^e.))] /kT)) (12) 

Eqs. ITT] and [T^ express the excess chemical potential as the sum of an average binding energy 
plus a fluctuation contribution involving deviations from the mean. 
Taking the average of Eqs. [TT] and [T2] and rearranging, we obtain 

1 kT 

I^Ts = {^es)o + ^ i{£es) " (^e.)o) + ^ (In (exp[(ee.- < £es >) / kT]) - lu {eiCg[- {s < £es >o)/kT])^) 

(13) 

The first term on the rhs of Eq. [T3] is q(f)np, where 0„p is the average net potential at the 
center of an uncharged solute. Then we can associate the bulk hydration free energy, nf^ 
above, with the sum of the non-polar contribution yU^^jy and all of the terms on the rhs of 
Eq. [13] except the first. The second term above involves the difference of the electrostatic 
interaction energy sampled with full coupling and with no electrostatic interactions; it has 
a substantial magnitude but contains no contribution from the surface potential. This term 
will contain a local contribution that is ion specific and a far-field contribution that does not 
depend on the ion charge or size and could be well-handled with continuum models. The 
final term on the rhs involves an average of fiuctuation contributions. Since previous work 
has shown that free energy shifts due to the surface potential can be handled at the mean- 
field level jol, this term thus includes little or no contribution from the surface potential. 



Further, it exhibits ion specificity due to relatively local interactions [47[. If the interaction 
energy distributions were exactly Gaussian, then the average of the fiuctuation terms would 
be zero (below). 

Considering now the first qcpnp contribution to the free energy, the net potential 0„p is 
obtained by sampling with no charge on the particle. If we assume the particle is situated 
at the center of a large water droplet, this potential will contain contributions from waters 
nearby the particle and from any asymmetric charge distribution at the distant water surface. 
Since all waters between these two domains will be randomly oriented, they contribute 
nothing to the net potential. Thus we can express the net potential as [sl 



■'np 



\ + <\>sp (14) 
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where 0/p is the local potential contribution from nearby waters. Based on the above discus- 
sion, both (pip and (psp contain large but canceling quadrupole contributions, leaving behind a 
net potential that reflects differences in water dipole orientations between the solute-solvent 
boundary and the liquid-vapor interface. 

Extensive previous classical point-charge simulation work has shown that, for these mod- 
els, the value of the local potential (pip is positive and of magnitude roughly 8 — 9 kcal/mol-e 
28|-|30| . The previous simulations have been performed in periodic boundaries, so the aver- 



age potential at the center of the uncharged particle is exactly (pip, with no (psp contribution 
from a distant interface. Addition of the (psp term (—13.8 kcal/mol-e for the SPC/E model) 
results in the (pnp value of —5.7 kcal/mol-e for the neutral iodi de p article in Ref. j^. For 
similar uncharged solute particles, Rajamani, Ghosh, and Garde |30| observed a variation in 
(pip of about 1 kcal/mol-e for modest-sized particles (all but the smallest Li"*" ion), displaying 
the small ion specificity to this quantity. Interestingly, Ref. 2^ shows that negative (pnp 
values persist up to large particle sizes of magnitude 15 A (since the observed magnitude of 
(pip is smaller than that for (psp), and the (pnp values are relatively size-independent beyond 
5 A. At very large sizes, (pnp should converge to zero for hard core solutes (as employed in 
Ref. ji^), but this convergence does not occur for molecular-sized solutes. 
Performing a second-order cumulant expansion of Eq. [T31 we obtain 

f^Ts ^ {ees)o + I iiees) " (^es)o) + ^ {(Ssi) " (^L) J (15) 

where fees is the deviation from the mean. Ref. [1^ shows that the difference of the fluctua- 
tion terms is very small for cations but of sizable magnitude for anions (with (feg^) > (^£L)o)- 



Outside a transition region centered at a small negative charge, Ref. [28| also suggests that 
linear response theory is reasonably accurate separately for cations and anions, but with 
different slopes. The cation/anion difference is illustrated by the result in Ref. 4^ that a 



simple Gaussian estimate of the free energy is modestly accurate for both electrostatically 
coupled and uncoupled sampling states for cations, but is only sensible for the coupled sam- 
pling state for anions. These differences arise from the penetration of the water protons 
closer to the anions. 

Alternatively, the second-order expansion of the electrostatic part of Eq. M yields 



Eq. [in] displays the equivalence to a Born model with the addition of the mean-field term 
q(pnp- The Born radius is related to the potential fluctuations at the uncharged solute center, 
and is thus temperature dependent. It has been shown previously that the distribution of 
total electrostatic interaction energies is not accurately Gaussian, and the deviations from 
Gaussian behavior are due to relatively local interactions between the ion and nearby waters 



48 , 



VI. ESTIMATING 



Starting from the model of Latimer, Pitzer, and Slansky (LPS) [43|], Ashbaugh and Astha- 



giri [4J] showed that, by shifting cation and anion radii by a constant distance (one distance 
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for cations and one for anions), and adding one other parameter that mimics the free energy 
for inserting the neutral particle, 



,»^(pa.) = 4(l-i)(_l^ + _l^)+2A (17) 

an extremely good linear fit is observed between Born and experimental hydration free 
energies for a wide range of monatomic ion pairs. Fitting to the pair free energies largely 
cancels any contribution from the net potential (with the caveat of the slight size-dependence 
to (psp)- These results suggest that the observed deviation from Gaussian behavior can be 
handled with two radius shift parameters. The smaller computed radius shift for anions is 
due to the observed stronger hydration [28| of anions compared with cations. 

After shifting all the single-ion free energies in Ref. [3| due to the standard state cor- 
rection of -1.9 kcal/mol discussed above, the cation and anion data in the second column 
of their Table III are shifted by an average magnitude of 11.6 kcal/mol from experiment 
[25I (with small deviations of order 0.5 kcal/mol). The shifts are quite uniform along the 
entire cation and anion series. Here it is suggested that this shift implies a (pnp value of -11.6 
kcal/mol-e (—0.5 V), and we consider this an experimental value. The fact that this result 
for the net potential (pnp is so close to the calculated value of (psp for point-charge classical 
models jl, [g, 0] has led to some confusion in the literature regarding these two quantities 
{(psp and (pnp)] the quantum mechanical results presented below highlight the large difference 
between them. 



Ref. [4J| also noted that free energies calculated from simulation (in periodic boundaries) 
are shifted from the LPS estimates in one direction for cations and the other for cations; 
these shifts are due to 4>ip (due to simulating in periodic boundaries), and not to Also, 
when the (pip shift is taken into account, the predicted bulk hydration free energies (here 
closely related to nf^ above) are nearly identical for large cations and anions. This result 
shows that previous simple models used to predict the electrochemical surface potential, or 
4>np, by subtracting a Born bulk estimate from the real hydration free energy of a large cation 



(with a result of —0.3 V) hold some validity [19|. 



VII. PARTITIONING THE FREE ENERGY 

Dividing up free energy contributions due to different regions of space, or different inter- 
actions, has a long history. In recent years, two such divisions have been derived 
that start from exact statistical mechanical theory. An exact formulation is helpful, since 
then the accuracy of subsequent approximations can be assessed. 

The Quasichemical Theory (QCT) [i^ involves a spatial decomposition of the free energy 
into inner-shell and outer-shell components. This partitioning is effected by manipulations 
involving repulsive particles. Here we will call the potential of interaction between the 
repulsive particle and the solvent M (model potential); the original QCT utilized a hard 
sphere particle for M, but a 'soft-cutoff' version has been derived [53(1 that allows for analysis 
using repulsive particles with continuous potentials. As in Ref. jSJ], the PDT can be re- 
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written as 



/x^^ = A;Tln(exp(-M/A;T)) - kT In (exp(-M/A;T))o - A;Tln (exp(-e/A;T)) 



M 



where e is the full solute-solvent interaction energy, and the sampling in the last term is 
conducted with the repulsive M-particle included. The first term on the right is an inner- 
shell contribution, the second term is a packing contribution that is the free energy to insert 
the M-particle, and the last term is the outer-shell, long-ranged contribution that involves 
solute-solvent interactions with the solvent pushed away from the solvent a distance specified 
by M. Physically, the inner-shell term reflects local (chemical) interactions of the solute with 
the solvent, and is minus the work required to push the nearby solvent out to a length scale 
specified by the M potential. The packing term is the free energy of cavity formation, and the 
long-ranged term includes all other interactions conditioned on the lack of solvent molecules 
in the inner shell. 

More recently, the Local Molecular Field Theory (LMFT) approach of Weeks and cowork- 
ers |56| has been adapted to calculations of free energies [48|. In terms of electrostatic 
interactions, the LMFT approach is to partition the interactions following the Ewald pre- 
scription: 

1 erfc(?7r) erf(?7r) 

where rj~^ specifies the length scale for the partitioning. In our studies of ion hydration j^. 



47|,|48l], that length scale has been chosen to mainly involve the first hydration shell, so r]~^ 
4 — 5 A. The first term yields a local electrostatic contribution, while the second includes all 
distant (far-field) interactions. In the Ewald method, the second term is transformed to a 
k-space representation. 

Using this LMFT partitioning, and first partitioning out the van der Waals (vdW) con- 
tribution, the free energy can then be written exactly as 

= l^Tdw + kT In (exp(£e.,;oc/A;T))^^^ + kT In (exp(£e.,/ar/A;T)) (20) 

Note that the sampling in the two electrostatic terms is conducted on different potentials, 
the local potential for the local term, and the fully coupled potential for the far-field term. 
Eq. [20] is equivalent to 

/^'" = f^7dw + kT\n {eM-ees,ioclkT)),dw + ^Tln {exp{-e,,jar/kT))^^^ (21) 
It was found in our models of ion hydration in bulk water (periodic boundaries, free of 



interfaces) and in the slab geometry [9|, |48|| , that the far- field electrostatic contribution is 
Gaussian distributed to high accuracy. 

Using the above LMFT partitioning, Ref. [9| showed that the free energy shift for an 
ion moving from a bulk periodic boundaries situation to the slab geometry is exactly the 
same as q4>sp, where (psp is determined from the integral of the charge density in Eq. [2J This 
is a direct confirmation that Eq. M is appropriate, albeit with large canceling contributions 
from the intrinsic free energy and the surface potential. This approach of adding a surface 
potential contribution to the intrinsic hydration free energy obtained in periodic boundary 
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simulations has been used widely recently in comparing single-ion (real) free energies to 



experiment during force-field development [g, |40|, |57 



The above discussion of energetic partitioning also provides a direct indication of why 



■'sp 



is not measurable thermodynamically: while can be measured [25|], separation into the 
intrinsic hydration free energy and surface potential parts involves an extra-thermodynamic 
partitioning of the electrostatic portion of the free energy. This point was noted by Guggen- 
heim 111 and Randies and Schiffrin 20 



VIII. THE CLUSTER-PAIR APPROXIMATION AND THE HYDRATION FREE 
ENERGY OF THE PROTON 

In 1998, Tissandier et al. experimentally determined the absolute hydration free energy 
of the proton j25|. They reported the value = —265.9 kcal/mol, obtained by combining 
cluster hydration data with bulk hydration data in a clever way (the cluster-pair approxima- 
tion, CPA). Once an absolute hydration free energy is obtained for one ion, single-ion free 
energies for all other ions can be obtained from accurate thermodynamic data. 

Here we briefly summarize the CPA in order to analyze whether the surface potential 
appears in the experimental results, and if so, where. We use the recently derived elegant 



and exact QCT formula from Asthagiri et al. [51 



/i'^^ = -fcTlnKfp^ + kTlnpin) + /i^^V - ^f^w (22) 

where K^^ is the equilibrium constant for binding n waters to the ion X in the gas phase, pw 
is the bulk density of water, p{n) is the probability of observing n waters in the inner-shell 
specified by a chosen cutoff radius, /ixvF„ hydration free energy of the XWn cluster, 

and /i^y is the hydration free energy of water in water. This formula is correct for any chosen 
n value, so long as that n does not exceed a limit imposed by the size of the observation 
volume. A sensible choice for n would be the value closest to the mean number of waters in 
the chosen observation volume, n. Then making the approximation p{n) ~ 1 will not lead to 
serious errors in the free energy. If that approximation is made, then the previously derived 
'primitive' QCT formula is obtained 5^ . 



We collect the first, second, and last terms together to obtain an effective free energy for 
inserting the ion into an n-water cluster: 

/i^ = -kT In [K^'^p-^lp{n)\ - n/i^ (23) 

leading to 

/i^^ = + /ixV (24) 

The formula allows us to connect exact statistical mechanics to the CPA. 

The CPA starts by defining conventional hydration free energies as deviations from proton 
values. Here we consider monovalent ions and label cations as P and anions as N . Then 

'^^ = /i^^ - (25) 
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and 

ex,con , ex : , ex l'r,Ci\ 

/iTV = f^N + /^H+ (26) 

These conventional free energies are available from thermodynamically measured free energy 
tabulations. An important point to note is that the conventional ion hydration free energies 
then include no contribution from the surface potential. 

Similarly, the conventional free energies for hydrating the clusters can be defined as 

I^PWn ^ /^pV„ ~ /^H+ (27) 

and 

f^NW,^ = f^NW^ + /^H+ (28) 

Using the above results, the following formula can be derived: 

ex,con ex,con] _ 1 \ ex ex 1 , ex (r)Q\ 

H'NWn I^PWn \ — 2 l^^NWn P PWn\ /^H+ l^-'J 



1 

2 

which is the same as 



1 r 
2 



ex,con ex,con ( ex ex X] _ ^ \ ex ex 1 , ex (nr\\ 

f^N H'P \h''N,n P'P,nJ\ — r) [H'NW„ H'PWnl ~'~ K'^^J 



The quantities on the left can be measured thermodynamically and in cluster experiments. 
The first term on the right clearly goes to as the number of waters increases. Following the 
discussion of Kelly, Cramer, and Truhlar 5^, data for a wide range of ions can be plotted in 



,1 r ■ r ex,con ex,con-[ /o i ■ r ex^con ex^coni /o i~i i x j. i- 

the form: x-axis as [/i^ ~ l-^p J/2 and y-axis as [/^atvi/^ — l^pwn J/2- ■'"o^ cluster studies 
(involving multiple ions) with sizes up to n = 6 water molecules, what is observed is that 
these curves all cross at a common point on the y-axis (see Fig. 1 in Ref. [E^). The slope 
of the curves decreases in magnitude with increasing n: for large n the curves should be fiat 
at the value of the proton hydration free energy, since the first term on the rhs approaches 
for large n. Assuming that the crossing point location holds as n increases, then the full 
Ihs value is a constant for all n at the crossing point. This suggests that the value of the 
first term on the rhs, for this hypothetical ion pair, is for all n, leading to an accurate 
free energy estimate even for small cluster sizes. The result of this analysis is not simply 
a cluster free energy; rather the cluster data is used to tease out an estimate of the bulk 
single-ion value. The only 'approximation' appears to be the assumption of continuation of 
the crossing behavior up to larger n. 

Given the above discussion, we can ask whether the surface potential is included in the 
reported free energy of the proton. Asthagiri, Pratt, and Ashbaugh [g^] have previously 
concluded that is indeed included in the reported value of —265.9 kcal/mol [1^. As 
discussed above, the conventional ion hydration free energy values contain no contribution 
from the surface potential. Then we can see from the left side of Eq. [30] that, in the limit 
of large n, the surface potential contribution is carried in the difference of the cluster free 
energies (yU^^ — Below, ion hydration free energies in clusters with n = 105 and 

n = 511 waters will be examined. These clusters display surface potential values close to 
the result computed in the slab geometry; the computed hydration free energy of an I~ ion 
agrees quite well with the value reported by Tissandier et al. j2^. Thus the results suggest 
both that the experimental proton value includes a contribution from the bulk water surface 
potential and that the CPA assumption holds for larger clusters. 
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IX. MEASURED ENTHALPY AND ENTROPY CHANGES FOR ANIONS NEAR 
THE INTERFACE 



From the basic PDT forms of the excess chemical potential, Eqs. |9] and [101 the hydration 
enthalpy and entropy can be derived from the appropriate temperature derivatives. Here we 
will ignore pressure-volume effects. The enthalpy is then 

/i- = {s) + UsR (31) 

where 

UsB^ = m^^x - {Un)n (32) 

is the difference in the average energy of the solvent with the ion (labelled X) present and 
with the ion absent (solvent reorganization energy). The resulting entropy is 

s- = ^-Hn(exp[(e-(5))]/A:T)) (33) 

These results show that the influence of the surface potential from a distant interface is 
carried predominantly in the enthalpy, not the entropy, since this effect will show up mainly in 
the mean binding energy (e), with a much smaller contribution from the entropic fluctuation 
term and no contribution from the solvent reorganization term. Ben-Amotz, Raineri, and 



Stell [61[ showed that the fluctuation contribution to the entropy is always negative. As 
noted above, it has been shown ^ that the free energy shift for an ion moving from the bulk 
to the slab geometries can be accurately modeled with the first-order, mean-field (enthalpic) 
term. This further supports the notion that the surface potential impact on the entropy is 
small. 

Also, the development shows that, while the solvent reorganization term is required for 
estimation of the separate enthalpy and entropy, it cancels out exactly when compiling the 
free energy. The cancellation is an example of enthalpy-entropy compensation. This has 
been noted previously by Ben-Amotz, Raineri, and Stell j6l|. Since the free energy profile 
determines the driving forces for ions near the interface, care should be taken in associating 
the term 'driving force' with separate enthalpy and entropy contributions. Only those en- 
thalpy and entropy components that can lead to a net contribution to the free energy can 
result in a driving force. It can also be noted that the measured very small magnitudes of 
the bulk entropies of ion hydration (they are similar to those for the isoelectronic rare gases) 
result from a near cancellation of a (positive) solvent reorganization term and a (negative) 
fluctuation term. 

Recent experiments have measured the surface adsorption enthalpy and entropy changes 
for a large chaotropic anion, SCN~ approaching the water surface from the bulk [s^]. The 
experiments indicate negative values for both quantities. The enthalpic term (-2.8 kcal/mol) 
is larger in magnitude than the entropic contribution TAs^^ (-1-2 kcal/mol), leading to a 
net surface adsorption free energy of -1.6 kcal/mol. Several simulation studies have explored 



the enthalpic and entropic contributions [62|, |63|. Yagasaki, Saito, and Ohmine [62| plotted 



the interaction parts of Eqs. ED and [23] (that is the Usr term was omitted), and observed 
positive values for the enthalpy and entropy changes. Caleman et al. [63| and Often et al. 
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33l | , on the other hand, computed the full enthalpy and entropy profiles indicated by Eqs. ET] 
and [33 and observed negative values for both quantities, in agreement with experiment. As 
indicated above, changes in the water-water interaction term as the ion approaches the 
surface contribute no net driving force, but lead to the observed negative values for Ah^^ 
and As*^^. Those separate quantities can yield insights into the ion hydration environment. 
It is suggested below that the net potential (pnp provides one part of the net driving force 
that tends to enhance anion density near the surface relative to cations. 

X. COMPUTATIONAL METHODS 

The calculations reported here include classical molecular dynamics (MD) simulations 
and quantum chemical calculations on cluster configurations extracted from the classical 
simulations. The classical MD simulations were performed with the Tinker code [g^], and 



the quantum calculations were performed with the ORCA code |65l . |66| , both freely available 



The classical simulations employed the SPC/E water model and used ion parameters from 



Refs. 57, [63] (parameter set 2 from Ref. (63], with a = 4.09 A and e = 0.1912 kcal/mol for 
the I~ ion and a = 2.70 A and e = 0.1554 kcal/mol for the Na"*" ion). Clusters of n = 105 
and n = 511 water molecules were modeled. The ion (or neutral particle) was restrained 
by a strong harmonic force to a particular distance from the center of mass of the water 
cluster. An external half-harmonic potential was employed to ensure that any evaporating 
water molecules were reflected back into the cluster; the starting point for the half-harmonic 
potential was chosen several A from the point where the radial density profile decayed to 
zero, however, so there was little or no interaction of the water molecules with the bounding 
potential. A time step of 2 fs was used, and all electrostatic interactions were included. The 
clusters were equilibrated for at least 500 ps, followed by production runs of at least 2 ns. 

Initial quantum calculations were performed on the n = 105 water cluster. Configurations 
were taken from the MD simulations that employed the SPC/E water model [3|. The calcu- 
lations were performed using Density Functional Theory (DFT) with the BLYP functional 
and the cc-pVDZ basis set, and included all of the electrons. This theory/basis combination 
was chosen to allow for calculations on a large set of 1000 configurations. The accuracy 
of the calculations in relation to higher-level methods will be discussed in the results sec- 
tion. FoUowup calculations aimed at probing the variation of the net potential with cluster 
size {n = 242), theory level (B3LYP and Hartree-Fock), and basis set (Aug-cc-pVDZ) were 
subsequently performed. The foUowup calculations were performed on a smaller subset of 
configurations (typically 109), so the results contain somewhat larger statistical fluctuations 
in the net potential, and those fluctuations are reported below. 

XI. CALCULATING THE SURFACE POTENTIAL 

As discussed above, the electrostatic profile 0(z) across the water liquid-vapor interface 
can be obtained from Eq. [H or equivalently from [7(] 

(P{z) - <f){v) = - r E,{z')dz' (34) 

J V 



14 



V 

e 

o 




50 100 150 200 

Time (ps) 



250 



Figure 1. The average surface potential as a function of time (ps) computed for the SPC/E n 
cluster using the far-field Ewald potential. 



511 



where Eg{z) is the electric field obtained by integrating the average charge density from 
the vacuum up to z. Alternatively, (j){z) can be obtained by simply averaging the Poisson 
potential computed during ab initio molecular dynamics (AIMD) simulations 0, [3|. While 
these expressions are exact, averaging the charge density or potential at a particular location 
results in substantial statistical fluctuations, especially for the quantum case consisting of 
electrons and the contribution from the pseudo-potential. 

Here, an alternative approach is proposed based on the LMFT (Ewald-like) Coulomb 
partitioning. Consider an arbitrary point chosen near the center of mass of a large water 
cluster. Then assuming that point is far from the liquid-vapor interface, the average charge 
density sampled over long times is zero everywhere near the point. Thus the local charge 
contribution to the potential at the point is zero. This conclusion of course holds all the 
way up to near the interface where the average charge density starts to deviate from zero. 
Then the surface potential is simply obtained from the average of the far-field part of the 
electrostatic potential at the chosen point. This approach eliminates the large potential 
fluctuations that occur while averaging the Poisson potential in AIMD simulations. The 
far-field contribution to the potential is smooth, and rapid convergence is observed. 

Fig. 1 shows the convergence of the surface potential as a function of simulation time 
for the n = 511 SPC/E water cluster. The surface potential converges to —14.9 kcal/mol-e 
(—0.65 V), and is relatively well-converged by 50 ps. The corresponding surface potential for 
the n = 105 cluster is —14.7 kcal/mol-e. An advantage of working in the cluster geometry is 
that there is no ambiguity due to intricacies of the Ewald potential that integrates to zero 
over the simulation box in periodic boundaries [9|. The surface potential computed for the 
n = 511 cluster is more negative than the slab result of Ref. [9| by 1.1 kcal/mol-e (0.05 V). 
As shown below, the quantum mechanical result from the far-field potential for the n = 105 
cluster (using the DFT models discussed above) is +67.0 ± 0.7 kcal/mol-e (2.9 V) obtained 
from only 21 statistically independent configurations. The small standard deviation provides 
a clear illustration of both the smoothness of the far-field potential and the resulting rapid 
convergence. 
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Figure 2. The water oxygen radial density profiles (arbitrary units) for the n = 105 (a) and n = 511 
(b) clusters, computed with the ion uncharged (black) and charged (red). 



A similar strategy can be imagined for analysis of the surface potential obtained from 
an AIMD simulation. In periodic boundaries, the far-field term is handled in k-space, as 
was done in the classical simulations of Ref. [9|]. The Ewald parameter r] can be chosen so 
that the majority of interactions with distant nuclei occurs outside the range of the nonlocal 
part of the pseudo-potential, and thus the effective charge for each nucleus is simply the 
total charge minus the charge of the omitted core electrons. Thus any potential ambiguity 
involving the pseudo-potential is also removed. 



XII. RESULTS 

Results are first presented concerning classical point-charge models of the I~ ion in n = 105 
and n = 511 SPC/E water droplets. The LMFT length scale parameter was chosen as 
r]~^ = 5.0 A. The radial water density profiles for the uncharged and fully coupled ions are 
shown in Fig. 2. The solute particle was restrained to lie at the cluster center of mass. It 
is clear from these plots that turning on the charge leads to a substantial rearrangement 
of the local hydration structure; the density profile near the liquid-vapor interface is not 
appreciably altered by the transition, however. 

Initially, the n = 105 cluster was simulated at 250 K in order to prevent evaporation 
events from the surface. The freezing point for a cluster is substantially lower than the bulk 
value, and the freezing temperature of the SPC/E model is 214 K [g^. Visual inspection of 
animations of the cluster motion show clearly that it is in a liquid state at 250 K. Subsequent 
simulations were performed at 300 K, resulting in a nearly identical water density profile with 
little or no evaporation. As discussed below, the surface potential and free energy properties 
of the cluster did not change appreciably between 250 K and 300 K. The quantum calculations 
reported below were performed primarily on configurations taken from the 250 K n = 105 
simulation, with one set of calculations on the n = 242 cluster at 300 K. 

The classical results for the potential at the center of the uncharged ion (j)np are first 
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Table I. Data for the surface potential and free energy calculations. All energies are in kcal/mol 
and all potentials are in kcal/mol-e. All errors are on the order of 0.1-0.2 kcal/mol or kcal/mol-e, 
except the QM surface potential calculation, with an estimated error of 0.7 kcal/mol-e. The value 
for the vdW free energy contribution was taken from Ref. 0] (slab geometry). The vdW value at 
250 K was shifted downward by 0.6 kcal/mol since the hydration entropy of the vdW particle is 
roughly -12.7 cal/mol-K. The finite-size dielectric correction for the distant interactions with water 
assumed cluster radii of 9.0 A for the n = 105 cluster and 15.0 A for the n = 511 cluster, resulting 
in corrections of —18.2 kcal/mol and —10.9 kcal/mol, respectively. All results are for the particle 
constrained to the cluster center of mass, except for those labeled by 'surf, which are results for 
the neutral particle located 8.5 A from the cluster center for the n = 105 cluster. The bulk periodic 
boundary free energy (from Ref. 0]) is shifted by —q(j)sp as implied in Eq. O 



CM potential data 



Geometry 


(psp 




4'np 




n = 105 (250 K, CM) 


-14.1 


+6.6 


-7.5 




n = 105 (250 K, CM, surf) 






-3.6 




n = 105 (300 K, CM) 


-14.0 


+6.9 


-7.2 




n = 511 (300 K, CM) 


-14.0 


+7.5 


-6.5 




Bulk Slab (300 K, CM) 


-13.8 


+8.1 


-5.7 




No solute n = 105 (300 K, CM) 


-14.7 








No solute n = 511 (300 K, CM) 


-14.9 








QM potential data 


n = 105 (250 K, QM) 


+67.0 


-76.9 


-9.9 




n = 105 (250 K, QM, surf) 






-3.8 




Free energy data 


Geometry 




^^es,loc 


^^es,far 




n = 105 (250 K, CM) 


3.9 


-30.4 


-14.6 


-59.3 


n = 105 (300 K, CM) 


4.5 


-31.5 


-14.5 


-59.7 


n = 511 (300 K, CM) 


4.5 


-31.9 


-20.3 


-58.6 


Bulk slab (300 K, CM) 








-58.3 


Bulk PBC (300 K, CM) 








-71.5 



presented (Table 1). The results show that the classical surface potentials calculated in 
the cluster and slab geometries are consistent (with a shift to more negative values by 1.1 
kcal/mol-e or 0.05 V), even for a cluster as small as n = 105, a case where a flat water radial 
density profile is not obtained near the surface. The cluster results with the uncharged solute 
included are shifted slightly to less negative values relative to the data presented above for 
the case of no solute present. This small shift is likely caused by local rearrangements of 
waters in the first hydration shell; the far-field potential smoothly grows in magnitude with 
increasing radius and thus a small contribution from that local region can lead to a small 
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shift. The local potential (jyip is positive and increases somewhat with increasing cluster size. 
Nevertheless, it is clear that the cluster calculations capture the essential physics of the 
potential at the center of the neutral particle, with a slightly more negative net potential 
compared with the previous bulk slab calculations. 

It can also be noted that, from the n = 105 simulations at the two temperatures, there is 
only a very weak temperature dependence to both and (pip (and the resulting (pnp), with 
most of the change coming from (pip. The results above show clearly that a sizable surface or 
net potential can exist, even with very small temperature derivatives of the potentials. The 
same conclusion was previously reached by Sokhan and Tildesley [5|. 

Utilizing the LMFT partitioning methods developed in Refs. [ol, 113, 48|, the computed 
hydration free energies are also shown in Table 1. The computed excess chemical potentials 
at the cluster center include one finite-size correction, which is a Born model correction that 
extrapolates to the large-cluster limit {fi^lei = ~(1 ~ l/e)9^/2-RD, where Rd is the droplet 
size). That dielectric correction includes no contribution from the surface potential (which 
is already in the cluster result), and is expected to be accurate for handling the distant 
interactions with water. The resulting free energies are in excellent agreement with the 
experimental value of —59.3 kcal/mol [25|. It is encouraging that insertion of the ion into a 
cluster as small as n = 105 can result in a quite accurate free energy estimate. Importantly, 
these results clearly include contributions from the surface potential, and the agreement with 
experiment thus suggests that the tabulation of Ref. 2^ does include that effect for all of 
the ion free energy values. We also note that, in agreement with previous results the 
far-field electrostatic contribution is accurately Gaussian in the cluster geometry. 

It is also interesting that, for the n = 105 cluster, the local electrostatic contribution to 
the free energy becomes less negative at lower temperature. This free energy change implies 
a positive local electrostatic contribution to the entropy, an effect previously observed in 
Ref. ji^]. There it was suggested that the positive value arises from a competition between 
hydrophobic and hydrophilic hydration for the large (chaotropic) anion. Consistent with 
Ref. |47|, the results also provide further confirmation that the fluctuation contribution to 
the entropy is relatively local. From the results for the total free energy change between 
the two temperatures, a rough estimate of +1.5 cal/mol-K results for the hydration entropy 
(including the negative non-polar contribution). This very small, positive result compares 
well with the value obtained experimentally (0.2 cal/mol-K) in Ref. [1^. Using different 



ion models, Ref. [47[ obtained a result of —0.3 cal/mol-K for the 1~ ion, and the results 
were compared to the tabulation in Ref. jl^] , exhibiting some deviation. Since the measured 
results in Ref. jl^ are likely of higher accuracy, we can consider the agreement encourag- 
ing. These results illustrate the strong ion specificity in hydration entropies, since small, 
kosmotropic ions display negative entropies of large magnitude [i^l . 

For comparison, the hydration free energy of the Na"*" ion in the n = 511 cluster was also 
computed, with a result of —105.8 kcal/mol. This result is shifted from the experimental 
value of —103.2 kcal/mol [1^. The resulting computed Nal total is —164.4 kcal/mol, com- 
pared with the unambiguous experimental total of —162.5 kcal/mol j25|. It can be noted 
that the force fields developed in Refs. jH^] and [s^] used a water surface potential value of 
— 12.1 kcal/mol-e, which is somewhat smaller than the computed value of —13.8 kcal/mol-e 



for the SPC/E model [9|]. The force fields were fit to ion hydration free energies and entropies 
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in reference to the experimental results of Ref. [25|. Thus, the parametrized force fields may 
reflect that difference. Also, calculations in the cluster geometry may differ mildly from bulk 
or slab geometries (see discussion below). If both I~ and Na"*" results are shifted upward 
by 0.9 kcal/mol to reflect the pair discrepancy, final values of —57.7 and —104.9 kcal/mol 
are obtained, respectively. Compared with experiment these values are shifted up by 1.6 
kcal/mol and down by 1.7 kcal/mol, respectively. These symmetric shifts likely are a direct 
reflection of the slightly larger surface potential magnitude observed for the clusters vs. the 
slab geometry. 

In order to obtain a rough estimate of the influence of the various terms in Eq. [T3] on the 
driving forces for ions near the water surface, each electrostatic term was computed for the 
= 105 cluster at 300 K for the I~ situated both at the cluster center and at a distance of 
8.5 A from the center (near the dividing surface). The nonpolar contribution was estimated 
from the slab data in Ref. j^. We consider changes in four terms: the nonpolar contribution, 
the net potential term, the binding energy difference term, and the fluctuation term. The 
resulting free energy changes are -4.2, -4.2, +10.2, and +0.9 kcal/mol, respectively, leading to 
a net repulsion of +2.7 kcal/mol. (Ref. [67|, Fig. 10 iodide set 2, shows that the free energy 
profile, using these classical models, begins its rise to positive values inside the dividing 
surface, and the radius chosen here was at the dividing surface for the cluster, leading to 
the positive free energy change.) The nonpolar and net potential portions both contribute 
a driving force towards the surface, while the binding energy term yields a large repulsion. 
This repulsion is due both to the loss of part of the hydration shell locally, and reduction in 
the distant interactions as well. The fluctuation term yields only a small repulsion. These 
calculations were redone in the slab geometry, with the following results (listed in the same 
order): -4.2, -2.6, +7.2, and +0.5 kcal/mol (the near-surface position for these calculations 
was chosen 0.5 A inside the dividing surface). Again, the change in the average fluctuation 
term is small. 

Finally, quantum mechanical (DFT) results are presented for the distribution of the net 
potential, and the results are compared with the classical distribution. Configurations 
(1000 in total) were drawn from a 5 ns classical simulation of the n = 105 cluster at 250 
K, as was done in Ref. j^. The quantum net potential at the center of the neutral particle 
was computed for each configuration using the ORCA code (DFT-BLYP, cc-pVDZ basis). 
The distributions were also computed with the particle restricted to a radius of 8.5 A from 
the center of mass to examine the electrostatic environment near the cluster surface. Fig. 3 
displays the results. 

It should be noted here that the classical particle was removed from the water cluster 
for computation of the net potential. Thus the computed potential is that for an ion- 
sized cavity, and is not exactly the same as {e)o in Eq. [121 which is the interaction energy 
of the ion with the surrounding waters sampled with no charge on the ion. The simple 
rearrangement into the form {e)cav + ((^)o ~" {£)cav) then shows that the correct evaluation of 
(e)o involves the cavity term plus the change in the local part of the potential (the far- field 
parts canceling) incurred by turning the point charge into the ion. Thus {s)cav calculated 
here is a first estimate of the net potential (pnp- The physical assumptions are that the ion 
charge distribution is spherical, there is no overlap with the water electron density, and only 
electrostatic interactions are included. In reality, the ion charge density will be distorted due 
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Figure 3. Net potential distributions at the center computed with classical and quantum models 
for the n = 105 cluster. The black curves are for the cavity located at the cluster center, while the 
red curves are for the cavity near the cluster surface. The full curves are the quantum (DFT-BLYP, 
cc-pVDZ basis) results and the dashed curves are the classical SPC/E results. 



to the nearby waters and there will be some electron density overlap. This ambiguity does 
not arise in the classical calculation. 

For the particle situated both at the cluster center and near the surface, the distributions 
of the cavity potential are relatively consistent between the quantum and classical SPC/E 
results. (A recent study suggests a classical model dependence; the TIP5P model results 
in a near-zero average net potential [lo|.) As shown in Table 1, the average net potential 
for the particle at the cluster center is shifted downward by 2.4 kcal/mol-e for the quantum 
calculation compared with the classical SPC/E result (to a value of -9.9 ±0.2 kcal/mol-e). 
This result suggests that the negative net potential is indeed maintained when modeling 
water at the quantum mechanical level. The classical and quantum net potentials computed 
at the cluster surface are very close to each other, with the average quantum result for the 
net potential shifted downward by only 0.2 kcal/mol-e. It can be noted that, while the 
quantum and classical distributions are very similar near the cluster surface (and skewed 
in the positive direction), the classical distribution exhibits a shift and a broader spread 
than the quantum distribution at the droplet center. Any asymmetry in the distributions 
is due to local interactions, since the far- field contribution is Gaussian distributed [9|. Also, 
the range of the distributions is large, spanning 35 kcal/mol-e or 1.5 V; thus aggressive 
sampling (5 ns classical simulation here) may be an important factor. Collecting the results, 
the net potential changes moving from the cluster center to the surface are +3.9 kcal/mol-e 
(classical) and +6.1 (quantum). The results imply a net driving force pushing anions toward 
the surface and cations away. 

Followup calculations were performed to test the robustness of the computed cavity po- 
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tential. Using the same basis and cluster size, DFT-B3LYP and Hartree-Fock calculations 
resulted in values of — 10.0±0.5 kcal/mol-e and — 8.7±0.5 kcal/mol-e, respectively. A DFT- 
BLYP calculation with the Aug-cc-pVDZ basis set (that includes diffuse functions) yielded 
a value of —7.4 ±0.5 kcal/mol-e. Finally a DFT-BLYP calculation (using the cc-pVDZ basis 
set) for an n = 242 water cluster resulted in a cavity potential of —7.9 ± 0.7 kcal/mol-e. 
Finally, one calculation was performed at the DFT-BLYP level in which the classical particle 
was sampled with a purely repulsive (WCA) potential as opposed to the full Lennard- Jones 
potential. The computed cavity potential for that case is —9.2 ± 0.5 kcal/mol-e. 

Of course the computed quantum net potentials can depend on the sampling employed 
here, namely using the classical SPC/E model over long times. This choice of sampling was 
chosen in order to allow for robust sampling of a diverse set of independent configurations 
for subsequent analysis. There is ample evidence, however, from both experiment j69| and 



ab initio simulations [70|, lDJ that the SPC/E model accurately reproduces the orientational 



distributions of the water molecules near the liquid- vapor interface. Another observation is 
that, assuming the net potential results from the relative water dipole orientations at the bulk 
and solute-solvent interfaces, the SPC/E result of -7.5 kcal/mol-e (for the n = 105 cluster) 
produces the quantum mechanically computed result of -9.9 kcal/mol-e when the classical 
result is scaled by the ratio 3.1/2.35; this is the ratio of the water dipole moments estimated 
from ab initio simulations j72| and the SPC/E dipole moment (2.35 D). Presumably, one 
component of the difference between the quantum result computed here (-9.9 kcal/mol-e) and 
that estimated from the LPS analysis above (-11.6 kcal/mol-e) is the above-noted difference 
between the cavity potential and the (5)0 quantity discussed above. 

In a recent ab initio CP2K-DFT simulation study (BLYP-D and PBE-D functionals) ex- 
amining the cavity potential in the bulk slab geometry [13], average potentials very close 
to or slightly positive were obtained, differing from the results presented here. Differ- 
ences between these two sets of calculations include: cluster vs. slab geometries, all-electron 
vs. pseudopotential calculations, Gaussian vs. dual Gaussian/plane- wave basis sets, and the 
means of sampling (classical SPC/E over 5 ns time scale vs. AIMD over 20-30 ps). Presently, 
the origin of the difference in the net potentials is not known. One source of potential con- 
cern is the incorrect asymptotic form for the potential in the DFT-BLYP model (and other 
gradient-corrected functionals) [t^ ; since the net potential results from the molecular dipole 
contributions of two different interfaces, and the potential may depend sensitively on the tails 
of the electron distributions, this issue deserves further investigation. As a preliminary step, 
the Hartree-Fock results presented here possess the correct asymptotic behavior. A result 
of approximately zero for the net potential would be somewhat surprising, since this implies 
an exact cancellation of the molecular dipole contributions between the bulk surface and the 
boundary between the small solute and water. Leung observed a close correspondence 
between the SPC/E and DFT results for the dipole contribution to the surface potential, 
but both models yielded a very small magnitude value due to the small system size (slab 
geometry). Naively, one might expect a higher dipole density near the uncharged solute due 
to the increased water density there. The classical (TIP4P-FQ) molecular dipole surface 
potential contribution profile in Ref. is positive, and if the dipole density is somewhat 
higher near the solute, then a net negative potential at the cavity center would be predicted, 
as we see in the SPC/E model. 
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Finally, results are presented for the quantum local and far-field contributions to the 
potential at the cavity center. To compute the far-field contribution, the charge density 
was obtained for 21 statistically independent configurations and written to cube files. The 
grid spacing for the electron density was taken as 0.08 au. On this grid, integration of 
the electron density led to slight deviations (on the order of 0.05 e) from a total charge of 
—1050. Thus all of the nuclear charges were scaled to yield a total positive nuclear charge 
that exactly balances the calculated electronic charge. Then, by comparison with individual 
net potentials computed with the ORCA code, the net potentials agreed to better than 1 
kcal/mol-e. It is expected that the deviations in the net potentials arise mainly from local 
contributions, so these results suggest the far-field values computed on the grid are of high 
accuracy. The resulting average far-field potential is -1-67.0 ibO. 7 kcal/mol-e (-1-2.9 V). Leung 
jsl has shown that the vast majority of this large value comes from the molecular quadrupole 
contribution. 

Two previous DFT studies of the water surface potential in the slab geometry have ob- 
tained values of -1-3.1 V [1] and -1-3.6 V [1]. The value reported here for the cluster geometry 
shows modest agreement with these estimates, but is slightly smaller in magnitude. The 
slight discrepancy likely arises due to the use of a modest-sized basis set in the present cal- 
culations. As discussed above, the potential jump across the interface involves a contribution 
from the trace of the molecular quadrupole moment. The quadrupole trace for the SPC/E 
model is -1-8.4 au (taking the center of mass as the origin) and the dipole moment is 2.35 D. 
The quantum mechanical trace for a single water molecule in the SPC/E geometry, using 
the BLYP functional and the cc-pVDZ basis is -38.1 au, while the dipole moment is 1.79 D 
(compared with the experimental value of 1.85 D). Using a higher-level aug-cc-pVTZ basis 
leads to a quadrupole trace of —42.3 au for the same geometry and a dipole moment of 1.74 
D. From these results it is clear that the large difference in (psp computed classically and 
quantum mechanically arises mainly due to the different signs and magnitudes of the water 
quadrupole moments in the point-charge and quantum models. A recent study has discussed 
the importance of the large quadrupole of the water molecule [t^], and alternative classical 
multipole and smeared-charge based water models can accurately reproduce the quantum 



quadrupole moment |75H77 



If the SPC/E structure is then minimized, and the aug-cc-pVTZ basis is used, the quan- 
tum quadrupole trace is —42.3 au and the dipole moment is 1.81 D. Calculations at the 
MP2 level for the relaxed geometry, using the same high-level basis set, result in virtually 
the same quadrupole trace and a dipole moment of 1.86 D. Thus we can see that the basis 
set used here underestimates the quadrupole trace by roughly 11%. Relaxing to the min- 
imum energy geometry does not substantially alter the quadrupole trace but does modify 
the dipole slightly. Of course in the condensed phase, water molecular polarization can lead 
to changes in both the dipole and quadrupole. If we simply scale the observed result of 2.9 
V by 1.1 (since the quadrupole contribution completely dominates the dipole contribution), 
however, the final result of 3.2 V is obtained, in decent agreement with the previous bulk 
results. The present cluster results are computed at the all-electron level and do not involve 
Ewald calculations of the potential. The agreement with the bulk results is therefore encour- 
aging and the results provide further confirmation of the large positive value of the surface 
potential when computed quantum mechanically. The present result also can be compared 
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with the value measured in electron holography experiments (+3.5 V) 

Taken together, the quantum <^„p result of —9.9 kcal/mol-e and the (^sp result of +67.0 
kcal/mol-e imply a (pip value of —76.9 kcal/mol-e. These results illustrate the large magnitude 
of the quadrupole contribution, as was previously shown for the bulk slab geometry [3|, and 
provide further circumstantial evidence for the complete cancellation of this contribution. 



XIII. SUMMARY DISCUSSION 



A principal aim of the present paper is to attempt to reconcile the now well established 
surface potential at the water liquid-vapor interface (3.1 to 3.6 V) with electrochemical 
estimates that fall in the range —0.4 to +0.3 V. The resolution comes through the rear- 
rangement of Eq. for the real hydration free energy into Eq. [71 The same analysis has 
been previously applied to classical studies to exhibit the cancellation of model-dependent 
quadrupole contributions that arise in the local and surface potential terms that sum to the 
net potential j^, 2^. The present paper has extended that analysis to quantum mechanical 
calculations. The quantum calculations are required for comparison to experimental surface 
potential measurements due to the very different classical and quantum water quadrupole 
moments. But we can conclude, based on the cancellation involved in Eq. [71 that the water 
molecular quadrupole contribution to the interface potentials does not play a role in ion 
hydration thermodynamics. 

While Eq. [S] is a correct expression of the electrochemical potential, it was shown above 
that the resulting intrinsic hydration free energies differ widely from any sensible bulk hydra- 
tion free energy estimates reported in the various tabulations. It is suggested here that the 
re-analysis of the Latimer, Pitzer, and Slansky (LPS) model ji^ by Ashbaugh and Asthagiri 
j3] provides the necessary connection to relate the disparate surface potential values arising 
from electron holography and electrochemical measurements. The single-ion LPS results ob- 
tained by fitting Eq. [TTjare a close approximation to the bulk hydration free energy defined 
above, 



If we compare the resulting bulk free energies to those measured in Ref. [25|, with a 
charge- dependent shift of magnitude 11.6 kcal/mol-e, excellent agreement is obtained for 
both the cation and anion series (to within less than 1 kcal/mol for each ion). We take 
this as an experimentally based estimation of 0„p = —11.6 kcal/mol-e (—0.50 V). Due to 
the consistency across the cation and anion series, we can thus consider the derived (pnp 
value to be mainly a property of water itself (at the bulk surface and near the spherical 
solute, with a small solute size dependence). Considering the relatively low level of the 
quantum mechanical calculations employed here, the cavity potential result obtained above, 
—9.9 kcal/mol-e (or —0.43 V), agrees modestly well with the LPS result. The results further 
suggest a i^l^ value for the proton of —254.3 kcal/mol to go with the real hydration free 
energy of —265.9 kcal/mol. This /i^^ value for the proton matches exactly the value given 



by Marcus (26| (when corrected for the standard state). (His value was determined using a 



method developed by Halliwell and Nyburg [78|; that method fits deviation of cation and 



anion free energies to a power series in the ion radius, starting with the third power, and 
extrapolating to large ions.) 
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A completely different approach to single-ion free energies was taken by Schmid, Miah, 
and Sapunov [27|. Their analysis rests on the assumption of equal aqueous entropies for 
the proton and hydroxide ions; the data in Ref. 25| result in a difference of less than 
3 cal/mol-K for these two ions, supporting that assumption. Following that assumption, 
their method uses only bulk pair hydration free energy data, and thus includes no interfacial 
potential contributions. Their result for the proton hydration free energy is —251.4 kcal/mol, 
shifted by 14.5 kcal/mol from Ref. [1^ (indicating 0„p = —14.5 kcal/mol-e). Alternatively 
the shift in the proton hydration enthalpy estimated from the TATB hypothesis yields a 
predicted net potential of —11.0 kcal/mol [79|. As can be seen, this relatively wide range 
of approaches yields a consistent indication of a shift between real and bulk hydration free 
energies, indicating a negative net potential 0„p. Further work is needed, however, to go 
beyond the estimation of the cavity potential presented here and accurately determine the 
(^es)o quantity in Eq. [121 Several recent papers have addressed the issue of the absolute free 
energy of hydration of the proton [13, 59, [g^, 80 -S^]- It is hoped that the present paper 
will help clarify the various discussions in those works concerning the role of the surface 
potential. 



The (t>np result discussed above differs from the conclusion of Randies ll8| based on the 



temperature dependence of cell emfs (the measured temperature derivative is small in mag- 
nitude and negative). Interpretation of those experiments assumes a smooth extrapolation 
to the critical point, leading to very small positive estimates of the surface potential. As 
noted in the results above, a small temperature derivative to the surface potential does not 
necessarily imply a small value for the surface potential itself, however. This point was 
previously noted by Sokhan and Tildesley j^. Thus we consider the negative value for (pnp 
obtained here as more convincing since it allows for reconciliation of bulk thermodynamic 
and recently obtained cluster data. 

In terms of implications for ion specificity in hydration, the resulting fif^ values for the K"*" 
and F~, which are of comparable size, are —74.4 and —116.0 kcal/mol, respectively. These 
values indicate a distinct charge asymmetry and show that anions are more strongly hydrated 
for a given ion size. Another possible implication relates to the recent debate concerning 



the relative propensities for the proton and the hydroxide ions for the water surface |84|-|86 



The measured isoelectric point for the surface is in the pH range of 2-4 (85|; if an average 
potential of roughly half the (pnp values above is assumed in the layer near the surface, the 
predicted protonation shift for water would fall in this range. 

Another conclusion of the present paper is that the influence of the water surface potential 
on ion distributions is related to the above discussion of bulk hydration free energies, and 
is not a separate issue j^, 3^. In fact, a part of the driving force selecting for anions at the 
water surface arises from exactly the same q4>np h'ee energy term. Eq. [12] makes the necessary 
connection, showing that the negative net potential provides a net driving force that moves 
anions to the surface while repelling cations. In colloquial terms, the ions never experience 
the full surface potential (psp, but rather (pnp, near the interface. Previous classical simulation 
results showing ions movin g ba ck into the bulk when the charge is changed from negative to 
positive support this view [GSl . 87, 8^. In addition, if the net potential effect were ignored, 
the Na"*" ion would be preferred at the water surface over the I~ ion . The net potential (pnp 
is due to long-ranged interactions, arising from the difference of surface and solute-solvent 
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boundary dipole layers, but generates a driving force on ions only within less than 10 A from 
the dividing surface. Of course other more local factors are crucial too in determining the 
details of the ion distributions near the surface. Taken together, the analysis of the hydration 
free energy shifts and surface-potential-related driving forces provide a more unified picture 
of ions near the interface. These conclusions leave one outlying case, however, which is that 
the F~ ion is repelled from the surface [s^, even with the driving force arising from the 
surface potential. The repulsion arises from detailed local interactions that differ near the 
surface, and will be the subject of a future paper. 

To close, a few observations are made concerning the results in the present study. One is 
that the very simple classical point charge models (and more accurate polarizable models) 
can yield sensible physical results for ions near the water surface, although the results are of 
course not identical to the quantum models. The results suggest that the potential outside 
of the classical water models can be relatively accurate in an average sense, while the charge 
distribution inside the molecules is completely unphysical [2|. Recent results suggest, how- 
ever, a model dependence to the cavity potential between the SPC/E and TIP5P classical 
models (see discussion above). 

The free energy results obtained from the classical droplet models are surprisingly robust, 
and provide further indications that the classical models represent a large fraction of the 
important physics of ion hydration. Neglected are more accurate representations of dispersion 



interactions |90l.l91|. polarization [89|,|92j, and chemical effects such as charge transfer [931.194 
effects that may have important impacts. A suggestion, perhaps heretical, is that modeling 
free energies in the cluster geometry may provide a more physically grounded approach for 
ion force-field development. One reason is that, instead of multiple periodic-boundary finite- 
size corrections there is only one free energy correction for the interaction of the 
ion with distant waters. This correction is likely accurate so long as the droplet size is greater 
than roughly 10 A. Another possibility is a general intuition that the physical environment 
at the center of a droplet may be more 'natural' than that in periodic boundaries. Observed 
shifts of local electrostatic contributions passing from the bulk (no interfaces) to the slab 
geometry provide a_hint of slight structural rearrangements near the ion with the introduction 
of free interfaces 



The present paper has focused on single ions in water droplets. The discussion has 



suggested that the water surface potential is included in the measurements of Ref. [25 
that the net potential rationalizes observed shifts compared with bulk estimates, and that 
the net potential influences ion distributions near the surface. A previous paper by Zhou, 
Stell, and Friedman [93], carried further by Pratt js^], provides a very different picture for 
the contact potential between two (conducting) bulk phases, each of which contains some 
(perhaps small) density of ions. Those papers show that the contact potential (psp is actually 
determined by the differences of the intrinsic hydration free energies of the ions between 
the two phases (due to charge neutrality in the two phases). The derivation can also be 
performed starting from Eq. [71 leading to an expression for the net potential (pnp in terms of 
the bulk hydration free energies fil^ discussed above. The net potential is not strictly ion- 
independent, but shows only small variation with ion size. Thus, as opposed to the discussion 
above focused on the effect of the surface potential on the ion distributions, the bulk ion 
hydration thermodynamics determines the surface or net potential. Further discussion of 
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these important results is outside the scope of this paper, but the conclusions illustrate the 
extreme subtlety involved in ion hydration thermodynamics. 
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